Set up directory, install packages and load packages

Exercise 1: Create study area map using “ggplot”

All datasets are found in the “data”-folder.

Here, we use the package ‘rnaturalearth’ to download countries in a ESRI shapefile format. These are available for download here: https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/

Coordinate reference systems (CRS)
CRS provide a standardized way of describing locations to describe geographic data. The CRS that is chosen depends on when the data was collected, the purpose of the data, etc. It is necessary to transform vector and raster data to a common CRS when data with different CRS are combined.

Setting CRS to SWEREF99 (Sweden). CRS can be referenced by its:

  1. EPSG code, CRS(“+init=epsg:3006”) (see http://www.epsg-registry.org/ and http://spatialreference.org/) or by

  2. proj4spring; “+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs”

There are two general options:

  1. unprojected (a.k.a. Geographic): Latitude/Longitude for referencing location on the ellipsoid Earth (wgs84), and

  2. projected: Easting/Northing for referencing location on 2D representations of Earth (the creation of maps) e.g., SWEREF99.

More reading:
https://www.nceas.ucsb.edu/sites/default/files/2020-04/OverviewCoordinateReferenceSystems.pdf

Load vector data

How to download using “rnaturalearth” and save the shp locally

I will open the saved .shp for the exercise:

## load country boarders -----------


#' downladed and saved from rnaturalearth (above)
swe_no_sweref <-  input_vector %>% 

  # value = T: return vector containing the matching elements
  grep(pattern = "swe_no", value = T) %>% 

  # read ESRI Shapefile object 
  sf::st_read() %>% 

  #convert geometry object into an sf object
  st_as_sf() 
## Reading layer `swe_no' from data source 
##   `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/country/swe_no.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -362704.5 ymin: -6097519 xmax: 1114145 ymax: 8972413
## Projected CRS: SWEREF99 TM
# check coordinate system 
st_crs(swe_no_sweref)$proj4string
## [1] "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# Create a bounding box as an sf object for layers ----
b_box <- c(xmin = 350432, ymin = 7179284, 
           xmax = 857780, ymax =7444230)

# Now intersect the full world's coastline with the bounding box 
suppressWarnings({ 
sweden.c <- st_crop(swe_no_sweref, b_box)
plot(sweden.c)
}) 

class(sweden.c)
## [1] "sf"         "data.frame"
### water data ----
input_vector
## [1] "input_data/vector/animal_data/data_cleaned_PV.rds"                                      
## [2] "input_data/vector/animal_data/data_ready_rsf.rds"                                       
## [3] "input_data/vector/animal_data/reindeer_data.rds"                                        
## [4] "input_data/vector/country/swe_no.shp"                                                   
## [5] "input_data/vector/country/swe.shp"                                                      
## [6] "input_data/vector/HR/spatial_files_availability_PV_winter_2021_shorter_modified_new.shp"
## [7] "input_data/vector/ne_10m_lakes/ne_10m_lakes.shp"
# import lakes
lakes <-  input_vector %>%
  grep(pattern = "lakes", value = T) %>%
  sf::st_read() %>%
  st_as_sf() 
## Reading layer `ne_10m_lakes' from data source 
##   `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/ne_10m_lakes/ne_10m_lakes.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1355 features and 41 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -165.9656 ymin: -50.66967 xmax: 177.1544 ymax: 81.95521
## Geodetic CRS:  WGS 84
class(lakes)
## [1] "sf"         "data.frame"
st_crs(lakes)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
# Transform or convert coordinates of simple feature (sf)
lakes <- st_transform(lakes, crs = 3006) 
suppressWarnings({ 
  
# crop an sf object to a specific rectangle (here based on extent of "sweden.c"-shp)
# sf_use_s2(FALSE) # use if errors occur below
lakes <- st_crop(lakes, sweden.c) 

class(lakes)
plot(lakes)
}) 

Create points of interest (cities)

Create xy points with labels. Here creating sf-object and setting the crs.

# Create points of interest (coordinates from QGIS or google maps for demonstration)
places <- data.frame(ID= c("Arjeplog", "Arvidsjaur", "Sorsele"),
                     y = c(7328407.945, 7280856.153, 7270326.684),
                     x = c(630040.101, 693262.414, 617477.638)) %>%
  
  st_as_sf(coords = c("x", "y"),  # create sf object 
           
           crs=st_crs(sweden.c))  # set same crs as for sweden.c     
class(places)
## [1] "sf"         "data.frame"
 st_crs(places) # check projection 
## Coordinate Reference System:
##   User input: SWEREF99 TM 
##   wkt:
## PROJCRS["SWEREF99 TM",
##     BASEGEOGCRS["SWEREF99",
##         DATUM["SWEREF99",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4619]],
##     CONVERSION["SWEREF99 TM",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",15,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Topographic mapping (medium and small scale)."],
##         AREA["Sweden - onshore and offshore."],
##         BBOX[54.96,10.03,69.07,24.17]],
##     ID["EPSG",3006]]

Make plot using ggplot

library(ggspatial)
plot(sweden.c)
## Warning: plotting the first 9 out of 168 attributes; use max.plot = 168 to plot
## all

plot(lakes)
## Warning: plotting the first 9 out of 41 attributes; use max.plot = 41 to plot
## all

plot(places)

ggplot() +
  
    # plot sweden and lakes 
    geom_sf(aes(geometry = sweden.c$geometry)) + # note: geom_sf() for vector data 
  
    geom_sf(data = sweden.c,                     # plot sweden border vector
          aes(geometry = sweden.c$geometry),
          fill= "antiquewhite")+

      geom_sf(data = lakes,                      # plot lakes vector
              fill="lightblue") +

  # layer_spatial(data=elevation) +
  
  # add cities 
  geom_sf(data = places) +
  
  # add city names 
  geom_sf_text(data = places, aes(label = ID),
               vjust=-0.5)+
  
  # Add text on map
  annotate(geom = "text", 
             x = 484609, y = 7404230, 
             label = "Norway", 
             fontface = "italic", 
             color = "darkgray", 
           size = 6) +
  
    annotate(geom = "text", 
             x = 784609, y = 7404230, 
             label = "Sweden", 
             fontface = "italic", 
             color = "darkgray", 
           size = 6) +
  
  # Add north arrow
  annotation_north_arrow(location = "bl", 
                       which_north = "true", 
                         height = unit(1, "cm"),
  width = unit(1, "cm"),
                       pad_x = unit(0.5, "cm"), # horizontal align
                       pad_y = unit(8, "cm"), # vertical align
                       style = north_arrow_fancy_orienteering) +
  #add annotation scale
  annotation_scale(location = "bl", 
                     width_hint = 0.5) +
  
  
  # keep new projection for ggplot
  coord_sf(crs = st_crs(3006), 
  datum = sf::st_crs(3006),             
  expand = FALSE) +
  
  # add grid lines 
  theme(panel.grid.major = element_line(color = gray(.5), 
                                        linetype = "dashed", 
                                        linewidth  = 0.1), 
        panel.background = element_rect(fill = "lightblue")) 

### examplefor removing or adding titles 
# theme(axis.title.x = element_blank(),    # remove 
#       axis.title.y = element_blank()) + 
  # xlab("Longitude") + ylab("Latitude")   # add

ggsave("output_data/Fig1_study_area.png", 
       width = 8, 
       height = 6, 
       dpi = 400)

Optional: Create a polygon and add to map

Example 2 interactive maps using mapview

# some gps data from reindeer 
reindeer_data <- readRDS("input_data/vector/animal_data/reindeer_data.rds") %>%
   dplyr::select(Collar_ID,  id, year, x_, y_, t_, group) 

str(reindeer_data)
## 'data.frame':    5951 obs. of  7 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ x_       : num  733706 733653 733836 733828 733683 ...
##  $ y_       : num  7314782 7314737 7314807 7314851 7314738 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : chr  "Fed" "Fed" "Fed" "Fed" ...
### basic usage ===========================================================

### most basic call
mapview()
### with some data
# mapview(reindeer_data) 
str(reindeer_data)
## 'data.frame':    5951 obs. of  7 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ x_       : num  733706 733653 733836 733828 733683 ...
##  $ y_       : num  7314782 7314737 7314807 7314851 7314738 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : chr  "Fed" "Fed" "Fed" "Fed" ...
# convert to sf 
reindeer_data_sf <- reindeer_data %>%
  st_as_sf(coords = c("x_", "y_"), #data frame to sf 
           crs=3006)              # assign to known crs 



plot(reindeer_data_sf)

mapview(reindeer_data_sf) 
# mapview structure
m <- mapview(reindeer_data_sf) 
m
str(m, 3) # two slots = object (points) and map (geometries)
## Formal class 'mapview' [package "mapview"] with 2 slots
##   ..@ object:List of 1
##   .. ..$ :sfc_POINT of length 5951; first list element:  'XY' num [1:2] 20.1 65.9
##   ..@ map   :List of 8
##   .. ..$ x            :List of 4
##   .. ..$ width        : NULL
##   .. ..$ height       : NULL
##   .. ..$ sizingPolicy :List of 7
##   .. ..$ dependencies :List of 5
##   .. ..$ elementId    : NULL
##   .. ..$ preRenderHook:function (widget)  
##   .. ..$ jsHooks      :List of 1
##   .. ..- attr(*, "class")= chr [1:2] "leaflet" "htmlwidget"
##   .. ..- attr(*, "package")= chr "leaflet"
m@object
## [[1]]
## Geometry set for 5951 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 19.75501 ymin: 65.85719 xmax: 20.24619 ymax: 65.92977
## Geodetic CRS:  +proj=longlat +datum=WGS84 +no_defs
## First 5 geometries:
## POINT (20.12738 65.86829)
## POINT (20.12616 65.86794)
## POINT (20.13028 65.86842)
## POINT (20.13018 65.86882)
## POINT (20.12679 65.86792)
class(m@map)
## [1] "leaflet"    "htmlwidget"
# mapview
mapview(sweden.c) # polygons
# character to factor needed for grouping colors
reindeer_data$group <- as.factor(reindeer_data$group) 
unique(reindeer_data$Collar_ID)
##  [1] "FIT_421952" "FIT_421943" "FIT_421941" "FIT_421940" "FIT_421949"
##  [6] "FIT_421910" "FIT_421939" "FIT_421946" "FIT_421906" "FIT_421956"
## [11] "FIT_421944" "FIT_421937" "FIT_421948" "FIT_421908" "FIT_421935"
## [16] "FIT_421942" "FIT_421954" "FIT_421907" "FIT_421932" "FIT_421931"
## [21] "FIT_421934" "FIT_421945" "FIT_421938" "FIT_421947" "FIT_421936"
## [26] "FIT_421953" "FIT_421909" "FIT_421950" "FIT_421951" "FIT_421933"
library(mapview)

mapid <- reindeer_data %>%
   filter(Collar_ID =="FIT_421952" | Collar_ID =="FIT_421943"| Collar_ID =="FIT_421941")  %>%
  st_as_sf(coords = c("x_", "y_"), crs=3006)
str(mapid)
## Classes 'sf' and 'data.frame':   628 obs. of  6 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : Factor w/ 2 levels "Control","Fed": 2 2 2 2 2 2 2 2 2 2 ...
##  $ geometry :sfc_POINT of length 628; first list element:  'XY' num  733706 7314782
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:5] "Collar_ID" "id" "year" "t_" ...
### Add home range 
AOI <- grep("modified_new.shp", input_vector, value = T) %>%
  sf::read_sf() %>%
  sf::st_transform(crs = 3006)
plot(AOI)
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf

pal <-  mapviewPalette("mapviewVectorColors") # optional: mapviewTopoColors, mapviewRasterColors, mapviewVectorColors, mapviewSpectralColors

mapid %>%
  mapview(zcol="Collar_ID", 
          cex = "Collar_ID", 
          col.regions = pal(3))+
  mapview(
    AOI,
    alpha.regions = 0,   # no fill
    color = "black",     # border color
    lwd = 2              # border width
  ) # or simply add + AOI 
### view extents
viewExtent(mapid)
mapview(st_bbox(mapid))
viewExtent(mapid) + mapview(mapid, zcol = "Collar_ID")
#### show tracks
mapid2 <- reindeer_data %>%
  filter(Collar_ID %in% c("FIT_421952","FIT_421943","FIT_421941")) %>%
  mutate(t_ = as.POSIXct(t_, format = "%Y-%m-%d %H:%M:%S", tz = "UTC")) %>%
  arrange(Collar_ID, t_) %>%
  st_as_sf(coords = c("x_", "y_"), crs = 3006)
str(mapid2)
## Classes 'sf' and 'data.frame':   628 obs. of  6 variables:
##  $ Collar_ID: chr  "FIT_421941" "FIT_421941" "FIT_421941" "FIT_421941" ...
##  $ id       : chr  "rt_PV_20_024" "rt_PV_20_024" "rt_PV_20_024" "rt_PV_20_024" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ t_       : POSIXct, format: "2021-01-29 20:00:44" "2021-01-30 00:01:20" ...
##  $ group    : Factor w/ 2 levels "Control","Fed": 2 2 2 2 2 2 2 2 2 2 ...
##  $ geometry :sfc_POINT of length 628; first list element:  'XY' num  733461 7315029
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:5] "Collar_ID" "id" "year" "t_" ...
tracks <- mapid2 %>%
  group_by(Collar_ID) %>%
 summarise(do_union = FALSE) %>%     # keeps point order within groups
  st_cast("LINESTRING")


mapview(tracks, 
        zcol = "Collar_ID", 
        lwd = 3, 
        layer.name = "Tracks") +
  mapview(mapid, zcol = "Collar_ID", cex = 3, layer.name = "Fixes") 
### basemaps
library(leaflet)
# To specify which background map you would like to use
mapview(gadmCHE, map.types = c("Stamen.Toner", "NASAGIBS.ViirsEarthAtNight2012"))
# see https://leaflet-extras.github.io/leaflet-providers/preview/ for more options -- to see which background maps that leaflet provides! 

Exercise 3 (optional): Maps in R

Pre-define crs using ESPG-code or proj4string by:

# projection to use as proj4string 
crs_to_use <- "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs" 

Load background data

Now let’s start by loading the maps using the “raster”-package. Another option: “terra”-package using rast()-function)

## [1] "input_data/ignore/PV_2021_cc.tif"    
## [2] "input_data/ignore/PV_2021_houses.tif"
## [3] "input_data/ignore/PV_2021_lichen.tif"
## [4] "input_data/ignore/PV_2021_nmd.tif"   
## [5] "input_data/ignore/PV_2021_roads.tif" 
## [6] "input_data/ignore/PV_2021_TRI.tif"   
## [7] "input_data/raster/PV_2021_dem.tif"   
## [8] "input_data/raster/PV_2021_slope.tif"
## [1] "Clear_cuts_2021"            "prox10_houses_stakke"      
## [3] "lichenmap_stakke_original"  "nmd10_stakke_ungeneralized"
## [5] "prox10_stakke_resample"     "ruggedness10_QGIS"         
## [7] "PV_2021_dem"                "slope10m_QGIS"
## [1] "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs"

## Reading layer `spatial_files_availability_PV_winter_2021_shorter_modified_new' from data source `/Users/hirn0001/Library/CloudStorage/OneDrive-Sverigeslantbruksuniversitet/Undervisning/Teaching and courses/GS_VMAS Data handling and illustrations/2025/create-maps/input_data/vector/HR/spatial_files_availability_PV_winter_2021_shorter_modified_new.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 716285.4 ymin: 7313268 xmax: 739293.8 ymax: 7321582
## Projected CRS: SWEREF99 TM
## [1] "+proj=utm +zone=33 +ellps=GRS80 +units=m +no_defs"

Add animal data (GPS points)

Load prepared reindeer data

## 'data.frame':    5951 obs. of  7 variables:
##  $ Collar_ID: chr  "FIT_421952" "FIT_421952" "FIT_421952" "FIT_421952" ...
##  $ id       : chr  "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" "rt_PV_20_020" ...
##  $ year     : int  2021 2021 2021 2021 2021 2021 2021 2021 2021 2021 ...
##  $ x_       : num  733706 733653 733836 733828 733683 ...
##  $ y_       : num  7314782 7314737 7314807 7314851 7314738 ...
##  $ t_       : chr  "2021-01-28 08:01:13" "2021-01-28 12:00:50" "2021-01-28 16:00:41" "2021-01-28 20:00:44" ...
##  $ group    : chr  "Fed" "Fed" "Fed" "Fed" ...
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"

Explore and select colors

## [1] "#00A600" "#01A600" "#03A700" "#04A700" "#05A800" "#07A800"

Plot and save all together

## Warning in plot.sf(HR_area, col = NA, border = 1, add = T): ignoring all but
## the first attribute

Some optional adjustments: hillshade instead of elevation

Prepare hillshade background map

Calculate slope and aspect from elevation (dem)

Plot and save all together

## Warning in plot.sf(HR_area, col = NA, border = 1, add = T): ignoring all but
## the first attribute